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SUMMARY 

Low  Cycle  Fatigue  (LCF)  in  engine  components  involves  macroscopic 
cyclic  plastic  strains  (with  a  stress-strain  hysteresis  loop)  over  a  significant 
portion  of  the  failure  region.  Characterising  elasto-plastic  behaviour  in 
potential  failure  regions  is  a  necessary  step  in  estimating  LCF  life. 

The  equations  governing  elasto-plastic  behaviour  are  summarized, 
and  the  methods  of  implementing  them  in  Finite  Element  (FE)  stress  analysis 
programs  discussed. 

An  extension  of  the  PAFEC  program  to  include  mixed  isotropic- 
kinematic  hardening  is  outlined,  and  verified  by  examples  for  which 
alternative  FE  solutions  were  available. 

A  sample  application  has  been  made  to  holes  in  a  plate  with  biaxial 
stress  fields  similar  to  those  in  disc  webs,  and  the  results  compared  with  the 
Neuber  and  modified  Stoweli  rules  commonly  used  for  design  life  estimation; 
these  rules  tend  to  overestimate  the  strain  in  biaxial  stress  conditions,  leading 
to  conservative  life  estimates.  /!  — 
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NOTATION 


General 

A  area 

c  kinematic  hardening  coefficient 

E  Young’s  modulus 

/  yield  function 

F  force 

g,  h  section  dimensions 

G  shear  modulus 

H  plastic  slope 

K  stress  or  strain  concentration  factor 
M  moment 

V  shear  force 

Vol  volume 

r  radius  or  radial  coordinate 

a  arc  distance 

x  linear  coordinate 

u  displacement 

a  rotation  of  face 

f)  rotation  derivative 

t  strain 

7  shear  strain 

k  shape  coefficient  for  shear  deflection 

A  proportionality  constant 

(i  kinematic  hardening  constant 

1/  Poisson’s  ratio 

< i>  angular  coordinate 

p  isotropic  fraction  in  mixed  hardening 

a  stress 

w  total  rotation 

Tensors 

C  elastic  or  plastic  constitutive  tensor 
J  tensor  invariant 

s  deviatoric  stress  tensor 

i  translated  stress  deviator 

a  translation  of  stress  origin 

6  Kronecker  delta 

(  strain  tensor 

a  stress  tensor 

Matrices  or  Vectors 

B  strain  -  displacement  transformation  matrix 
D  stress  -  strain  constitutive  matrix 
F  load  vector 

u  displacement  vector 

£  strain  vector 

#  stress  vector 


'i 


Subscripts  or  Superscripts 

A  applied  loads 
E  elastic 
P  plastic 

y  current  yield  condition 
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1.  INTRODUCTION 

As  engine  power  and  speed  is  varied  in  an  aircraft  turbine  engine,  vital  components 
are  subject  to  changes  in  thermal  and  mechanical  loading.  Each  flight  includes  one  major 
start -stop  cycle,  plus  lesser  cycles  depending  on  the  manner  of  operation  of  the  engine. 
High  stress  levels  are  used  in  many  components  to  achieve  lightness  and  consequent 
high  specific  output,  and  cyclic  variation  of  stresses  in  the  components  leads  ultimately 
to  cracking  at  critical  stress  locations,  with  failure  by  the  process  of  low  cycle  fatigue 
(LCF). 

There  are  no  practical  or  economic  means  of  determining  by  inspection  when  the 
LCF  life  of  components  is  exhausted,  so  accurate  and  reliable  means  of  life  estimation 
are  needed.  Three  categories  of  data  are  required  for  LCF  life  prediction: 

(i)  The  loading  history,  preferably  obtained  from  operational  records  for  the  engine, 
alternatively  derived  from  simulation  or  estimates. 

(ii)  Modelling  and  analytical  methods  which  use  the  operational  loading  to  produce 
the  history  of  stress,  strain  (and  temperature  where  relevant)  at  critical  locations  in 
components  where  failure  may  originate. 

(iii)  A  damage  accumulation  criterion,  so  the  cumulative  effect  of  loading  cycles  can 
be  quantified  and  a  measure  of  fatigue  damage  provided.  This  may  be  related  either 
to  the  life  to  crack  initiation  (for  safe  life  estimates),  or  to  estimated  crack  propagation 
rate  (in  damage  tolerant  design  using  fracture  mechanics  methods). 

Thermodynamic  and  heat  flow  analyses  of  engine  operational  data  are  used  to  de- 
. ermine  generalised  loadings,  including  pressures,  metal  temperatures  and  rotor  speeds. 
Finite  tlement  (FE)  models  handle  the  complex  boundary  shapes  commonly  found  in 
practical  components,  so  FE  computer  programs  are  generally  used  for  determination 
of  detailed  temperature,  stress  and  strain  distributions. 

Since  LCF  by  definition  is  characterised  by  macroscopic  cyclic  plastic  strains  mani¬ 
fested  by  a  stress-strain  hysteresis  loop  over  a  significant  portion  of  the  failure  region,  a 
material  model  incorporating  cyclic  plastic  behaviour  is  required.  Predicted  stress  and 
strain  values  depend  on  plastic  properties  of  the  material.  Further,  when  the  component 
is  unloaded,  plastic  behaviour  leaves  residual  stresses  which  also  influence  fatigue  life. 

High  cycle  fatigue  life  predictions  are  based  on  elastic  stress  estimates  (Basquin  law), 
but  strain  amplitude  is  regarded  as  the  best  parameter  for  predicting  LCF  life.  Total 
strain  values  (elastic  plus  plastic)  are  used  in  the  Coffin-Manson  equation1,  adopted 
almost  universally  for  LCF  prediction.  Hence  characterising  elasto-plastic  behaviour  in 
potential  failure  regions  is  a  necessary  step  in  estimating  LCF  life. 

This  report  summarises  the  equations  governing  elasto-plastic  behaviour,  indicates 
how  they  are  implemented  in  FE  stress  analysis  programs,  evaluates  and  compares  FE 
packages  available  at  ARL,  describes  the  extension  of  the  PAFEC  program  at  ARL 
to  include  mixed  isotropic  kinematic  hardening,  and  gives  applications  of  FE  elasto- 
plastic  analysis  under  loading  resembling  engine  components,  comparing  results  with 
the  commonly  used  Neuber  and  Stowell  approximations. 

2.  CYCLIC  PLASTICITY  ANALYSIS 

In  addition  to  equilibrium  and  compatibility  (strain-displacement)  relations  used 
in  elastic  stress  analysis  (and  not  changed  in  form  if  strains  are  assumed  to  be  small), 
three  concepts  are  required  to  formulate  cyclic  plasticity  problems2  . 

(i)  Initial  Yield  Criterion 

A  number  of  yield  criteria3  •  *  have  been  formulated,  but  only  the  maximum  shear 
of  Tresca.  and  the  distortion  energy  or  octahedral  shear  of  von  Mises  are  used  for  normal 
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ductile  metals.  For  the  more  commonly  used  materials  the  von  Mises  criterion  agrees 
better  with  test  data,  and  is  mathematically  simpler  as  it  defines  a  single  continuous 
failure  surface.  This  surface  is  expressed  in  terms  of  a  function  of  the  second  invariant 
(J2  =  »ij  *ij )  of  the  deviatoric  stress  tensor  «,y  and  can  be  written: 


where  <r0  is  the  initial  yield  stress  in  uniaxial  tension 

»ij  =  a ij  -  f>ij  «tjh/3  ,  with  stress  tensor  tr.y  and  Kronecker  delta  6tj 
Ukk  =  <?n  +  022  +  using  the  tensor  convention  that  repeated  dummy 
subscripts  indicate  summation 

In  3-dimensional  principal  stress  space,  the  von  Mises  criterion  describes  a  cylindri¬ 
cal  surface,  of  radius  o0  =  (§s(J  «,/) »  with  axis  through  the  origin  and  equally  inclined  to 
the  three  principal  stress  axes  (Fig.  1).  For  plane  stress  (<73  =  0),  this  surface  intersects 
the  01,02  principal  stress  plane  in  an  ellipse. 

(ii)  Plastic  Flow  Law 

Total  increments  in  strain  are  assumed  to  be  divisible  into  elastic  and  plastic  frac¬ 
tions: 

dtij  =  dtfj  +  dt^ 

The  direction  of  the  plastic  strain  increment  tensor,  according  to  the  associated  flow 
rule,  is  normal  to  the  yield  surface  (requiring  maximum  plastic  work  on  deformation): 


=  dX 


9f 

doij 


where  dX  is  a  constant  of  proportionality.  This  is  the  normality  condition  and  for  a  von 
Mises  material,  is  equivalent  to  the  Prandtl-Reuss  equations  =  <IX  . 

(iii)  Strain  Hardening  Rule 

In  many  materials  yield  strength  increases  progressively  with  increasing  plastic 
strain.  While  the  von  Mises  yield  criterion  and  the  normality  flow  rule  for  plastic  strain 
are  generally  accepted  and  adequately  established  by  experiment,  a  similar  consensus 
does  not  exist  for  rules  used  to  describe  strain  hardening.  Simpler  rules  are  preferred, 
provided  their  behavioural  description  is  reasonable,  as  more  elaborate  hardening  models 
lead  to  considerable  complexity  in  finite  element  programs.  Four  main  hardening  rules 
have  been  used,  viz.  isotropic,  kinematic,  Mroz  or  multiple  surface  models,  and  sublayer 
or  subvolume  models. 


(a)  Isotropic  Hardening 

The  simplest  assumption  is  that  material  strengthens  tiniformly  with  increasing 
plastic  strain,  irrespective  of  strain  direction  forward  and  reversed  loading  have  equal 
effects  (Fig.  2).  This  implies  that  the  yield  surface  radius  expands  uniformly  in  stress 
space,  maintaining  shape,  orientation,  and  origin  of  axis  (Fig.  3).  Isotropic  hardening 
does  not  satisfactorily  model  reversed  or  cyclic  loading,  but  its  implementation  is  simple 
and  for  monotonic  loading  results  are  as  good  as  for  more  complex  models5 . 

(b)  Kinematic  Hardening 

When  loading  is  reversed,  most  metals  exhibit  a  reduced  yield  strength  in  the 
reverse  direction,  known  as  the  Bauschiuger  effect  (Fig.  2).  The  kinematic  hardening 
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model  formulated  by  Prager®  represents  this  effect  by  assuming  that  the  yield  surface 
translates  as  a  rigid  body  in  stress  space,  maintaining  its  size,  shape,  and  orientation 
(Fig.  4).  This  translation,  da,/  is  in  the  direction  of  the  plastic  strain  increment,  i.e. 
normal  to  the  yield  surface.  The  Prager  formulation  is  not  invariant  in  a  subspace 
of  reduced  dimensions  -  it  must  be  modified  for  a  two  dimensional  system.  Ziegler15 
formulated  an  alternative  definition  of  the  yield  surface  translation,  da,/  namely,  that  its 
direction  is  the  vector  connecting  the  centre  point  (or  origin)  of  the  current  yield  surface 
to  the  existing  stress  point.  In  most  problems,  the  difference  in  the  results  from  the 
Prager  and  Ziegler  formulations  is  negligible,  with  the  Prager  equations  generally  more 
convenient  to  use.  Both  satisfactorily  model  reversed  loading  with  elastic-perfectly 
plastic  behaviour,  or  when  the  degree  of  strain  hardening  is  limited,  and  hardening 
is  reasonably  linear.  Difficulties  arise  at  large  reversed  strains  when  the  stress-strain 
relation  for  the  material  is  highly  nonlinear,  because  there  is  no  satisfactory  rule  for 
relating  hardening  coefficient  to  cumulative  strain. 

(c)  Mixed  Isotropic-Kinematic  Hardening 

Mixed  hardening  formulations 7 ' 8  have  been  developed  combining  isotropic  and 
kinematic  hardening,  which  means  that  the  yield  surface  both  expands  and  translates. 
By  extending  the  scope  of  the  kinematic  model,  better  agreement  with  test  data  is  ob¬ 
tained  for  some  materials.  A  particularly  simple  and  readily  implemented  formulation 
of  mixed  hardening  has  been  developed  and  is  presented  below. 

(d)  Multi-Surface  Models 

Non-linear  hardening,  reduced  in  the  uniaxial  case  to  a  piecewise  linear  stress-strain 
curve,  is  represented  by  the  Mroz  and  subvolume  models.  The  Mroz  model 9  comprises 
a  series  of  (initially  concentric)  yield  surfaces,  each  related  to  a  particular  yield  stress 
corresponding  to  one  segment  of  the  curve.  When  plastic  strain  occurs,  the  surface 
translates  until  it  touches  the  next  bounding  surface,  which  is  translated  in  turn.  Contact 
between  the  surfaces  is  maintained  until  unloading  occurs  (Fig.  5). 

A  similar  piecewise  linear  modelling  is  achieved  by  the  sublayer  or  subvolume 
models10,11  .  These  postulate  that  the  material  comprises  a  number  of  subvolumes 
(physically  analogous  to  material  grain  structure).  Each  subvolume  has  elastic-perfectly 
plastic  properties  with  a  different  yield  stress.  All  are  subject  to  the  same  strain,  with 
their  individual  stresses  combining  to  give  total  stress. 

The  Mroz  and  subvolume  models  result  in  the  Masing  description 12  of  the  Bausch- 
inger  effect,  viz.  that  on  reversed  loading  the  shape  of  initial  loading  curve  is  maintained, 
magnified  by  a  factor  of  two.  This  is  a  satisfactory  representation  for  many  engineering 
materials.  More  elaborate  models13,1®’23  have  been  formulated  to  describe  fine  de¬ 
tails  of  behaviour,  but  generally  these  are  too  complex  to  implement  in  a  finite  element, 
computer  program,  or  too  restricted  in  the  materials  represented. 

3.  PLASTICITY  EQUATIONS 
Failure  Criterion 

A  general  form  of  the  von  Mises  yield  criterion,  incorporating  the  various  hardening 


rules,  can  be  expressed  as: 

/  =  (j**/)  “  (T*(,P)  =  0 

(1) 

where 

C; 

II 

C; 

i 

B 

v.; 

(2) 

5 

and  Sij  =  Oij  -  <5,y  cr*  t/3  as  previously  defined.  <r,(tp)  is  the  yield  stress  as  a  function 
of  the  plastic  strain.  The  translation  a,y  of  the  stress  origin  is  also  a  function  of  plastic 
strain,  depending  on  the  hardness  model  adopted,  and  determines  a,j ,  the  translated 
deviatoric  stress. 


Plastic  Flow  Equations 

Strain  Partitioning 

Normality  Condition 


dtH  =  dtf,  +  d*ij 


3  dX 
2  a. 


(3) 


(4) 


The  proportionality  constant  dX  is  found  by  taking  the  inner  product  of  equation  (4): 

2  ,  p  ,  p  _  3  d A2 

2  diij  dtij  ~  2  ff2  S‘J  S‘J 


so: 


d*  =  =  dt* 

For  a  von  Mises  material  dX  =  dtp  is  equal  to  the  strain  in  imiaxial  tension. 

Strain  Hardening 

Isotropic  Q,y  =  0 

da,  =  Hdtp 

where  H  is  the  plastic  slope  in  simple  tension,  so  a,  =  ay  ,  the  yield  stress  in  tension. 
Kinematic  da,  =0  so  a,  =  a0  ,  the  initial  yield  stress. 
daij  =  cdtp-  (Prager) 

daij  =  dfi(ai}  -  aij)  (Ziegler) 

Equivalence  with  uniaxial  tension  requires  that  the  constants  c  =  \H ,  or 
H  .  - 


(5) 


dfi  =  di  so  with  equations  (4)  and  (5): 


da,-.-  =  -±Hdtp 

a. 


or: 


(Prager) 

(G) 

ip  (Ziegler) 

(T) 

Mixed  Hardening 

The  hardening  effect  H  d(P  in  equations  for  isotropic  and  kinematic  hardening  can 
be  divided  into  fractions  p  defining  yield  surface  expansion  (isotropic),  and  (1  —  p ) 
yield  surface  translation  (kinematic).  Notionally  either  the  strain  hardening  rate  H 
or  the  plastic  strain  increment  dtp  can  be  so  partitioned,  and  this  makes  subsequent 
implementation  in  computer  code  particularly  straightforward. 


da,  =  pHdtp 

(8) 

so: 

a,  —  a 0  +  Pi17!)  ~  ffo) 

dan  =  ( I  —  p)Hdip  (Prager) 

v. 

(9) 
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(  :  dan  =  ^(1  -  p)Hdtp  (Ziegler)  , 

a. 

Constitutive  Equations 

Constitutive  equations  relating  increments  of  stress  and  elastic  strain,  using  equa¬ 
tion  (3),  can  be  written:  , 

dctij  =  Cfikldtfi  =  Cfjkl  (dckt  -  d(kl)  (10) 

where  the  elastic  stiffness  tensor  is: 

C?jkl  —  2 G(S{kSji  +  j  _  —S{ jSki) 


with  shear  modulus  G.  Poisson’s  ratio  v . 

Constitutive  equations  for  plastic  deformation  are  derived  by  noting  that  during 
plastic  flow  the  stress  point  stays  on  the  yield  surface,  so  that  df  —  0  in  equation  (1): 


df  =  IL-dan  +  j-dan  +  -Ida.  =  0 


doij  ’’  dan 


(11) 


Substituting  in  equation  (1)  and  using  (2),  (8)  and  the  Prager  result  of  (9): 

3  ~‘  {dan  -  dan)  ~  pHdtp  =  0 

2  <T , 

1  ^  {dan  -  (1  -  /»)  —  Hdep)  -  pHdtp  =  0 

2  ot  y  ct  f  ' 

with  equations  (3)  and  (10): 


(12) 


Since  in  =  0,  for  an  elastically  isotropic  material  only  the  “diagonal”  components 
i  =  k,  j  =  l  of  the  tensor  Cpkl  produce  non-zero  terms,  and  equation  (12)  can  be 
reduced  to: 

“  =  ^%gS  (13> 

Using  equations  (4)  and  (13)  with  (10),  with  dummy  subscripts  k,l  changed  to  m,n: 

3  8kl&mn  dCfj 


,  (.  ^  8kt9mndtmn  1 

d<r a  =  ('nki\  diki  -  2  ^(T+h/zc)  I 


dan  =  (Cfjkl  -  CPiki)dtki  (14) 

By  a  reduction  similar  to  equation  (13),  the  plastic  constitutive  tensor  becomes: 


_  3(7 

,]kl  ~  1  +  tf/3(7  a] 


(15) 
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4.  FINITE  ELEMENT  IMPLEMENTATION 


Finite  Element  Analysis 

For  a  continuum,  the  constitutive  equation  (14)  can  be  used  in  standard  finite 
element  equations  of  incremental  form.  Using  matrix  notation: 


/  BT(  D  -  Dp  )B  dVol  £u  =  ^ F.4 

(16) 

(K-Kp)4u  =  *Fa 

(17) 

where  K,KP  are  elastic  and  plastic  structural  stiffness  matrices 
£F.a  is  vector  of  applied  loads 

B  is  displacement  strain  transformation  matrix  so  St  —  B^u 
D,DP  are  material  elastic  and  plastic  constitutive  matrices, 
corresponding  to  tensors  C^kt  in  equation  (14). 

Equation  (17)  provides  an  incremental  solution  for  6u  using  the  tangential  stiffness 
matrix  (K  -  Kp).  Be  c  ause  Kp  is  a  function  of  stress,  the  solution  is  obtained  by  suc¬ 
cessive  Newton-Raphson  approximations,  and  the  stiffness  matrix  has  to  be  re-evaluated 
and  reduced  at  each  iteration. 

Repeated  matrix  reductions  can  be  avoided  by  re-arranging  equation  (16)  so  that 
only  the  constant  elastic  matrix  K  need  be  decomposed,  which  is  required  once  only 
during  the  solution: 


/  BtD  B  dVolSu  =  SFA  +  JBTDpBdVo/<5u  (18) 

K  =  tfF.A  +  SFp  (19) 

The  plastic  pseudo-forces  SFp  are  found  by  noting  that  when  equations  (10)  and  (14) 
are  written  in  matrix  form: 


=  (D  -  Dp)#<  =  D(<5t  -  Str)  (20) 

Hence  Dp6r  =  Dtfcp  and  since  St  =  B6u  the  RHS  integral  in  equation  (18)  becomes: 

SFp  =  f  BTD  dVolfitp  (21) 

Partitioning  i^u  =  Sup  +  Sup  in  equation  (19)  leads  to: 

K(*u£  +  Sup)  =  ^F,a  +  *Fp  (22) 

The  equation  K  flue  =  ^F.a  is  simply  a  scaled  elastic  solution  which  can  be  subtracted 
from  equation  (22)  giving: 

K  £u/>  =  AFP  (23) 

In  this  equation  />U/>  is  found  by  iterative  solution,  with  successive  approximations  to 
SFp  given  by  equation  (21). 

Computer  Programs 

The  commercial  FE  program  package  PAFEC  was  available  at  ARL  as  a  general 
purpose  stress  analysis,  heat  flow  and  dynamics  program  when  this  work  commenced. 
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Subsequently  access  to  NASTRAN  was  acquired  for  special  applications,  and  the  plas¬ 
ticity  capabilities  of  both  program  packages  are  now  briefly  compared. 

PAFEC  uses  the  development  just  presented;  a  flow  chart  of  the  PAFEC  plasticity 
routine  is  shown  in  Figure  6.  This  approach  is  simple,  but  at  times  it  results  in  conver¬ 
gence  difficulties.  NASTRAN  by  comparison  uses  the  tangential  stiffness  approach  of 
equation  (17),  though  a  flexible  strategy  is  employed  by  updating  the  matrix  only  after 
an  optimum  number  of  iterations.  Convergence  is  assured,  but  computational  effort  is 
greater. 

Regarding  types  of  elements  available  for  elasto-plastic  problems,  PAFEC  is  much 
more  versatile,  and  NASTRAN  is  highly  restricted  in  its  ability  to  model  continua  of 
complex  shape,  typical  of  aircraft  engine  components.  For  2-dimensional  models,  NAS¬ 
TRAN  has  only  straight-sided  linear  elements,  with  plastic  conditions  represented  only 
at  the  centroid.  In  three  dimensions  more  capability  is  offered  by  linear  brick  elements 
with  2x2x2  Gauss  point  integration,  and  equivalent  triangular  prism  elements.  By 
comparison  PAFEC  has,  as  well  as  linear  elements,  quadratic  and  cubic  isoparametric 
elements  in  two  and  three  dimensions,  with  2x2  or  2x2x2  Gauss  integration.  Both 
triangular  and  rectangular  (or  brick  and  prism)  shapes  are  available,  and  the  curved 
shapes  obtainable  with  isoparametric  elements  allow  better  modelling  of  practical  com¬ 
ponents. 

Hardening  rules  available  are  similar.  Both  programs  have  isotropic  and  kinematic 
models.  NASTRAN  also  offers  “mixed”  hardening,  but  this  is  restricted  to  the  special 
case  of  p  =  0.5  in  equations  (8)  and  (9),  which  implies  that  the  yield  stress  always 
remains  equal  to  the  initial  value  despite  reversal  of  loading.  The  writer  has  extended 
PAFEC  code  to  incorporate  the  more  general  case  of  mixed  hardening  described  by 
these  equations.  Nonlinear  hardening  is  available  with  the  isotropic  model,  but  with 
kinematic  hardening  PAFEC  expresses  the  yield  surface  translation  in  an  equation  of 
the  form  a,;  =  ,  which  requires  constant  H  (i.e.  linear  hardening)  when  integrating 

equation  (G). 

5.  COMPARISON  OF  SOLUTIONS 

As  a  preliminary  to  using  PAFEC  in  practical  analysis  of  cyclic  plasticity,  suit¬ 
able  verification  problems  were  sought.  The  only  reasonably  accessible  example  which 
included  linear  strain  hardening,  was  the  proving  ring  problem  in9  ;  an  alternative  so¬ 
lution  to  this  problem  showing  slightly  different  results  is  also  found  in15.  The  same 
problem  with  mixed  hardening  is  treated  in7  .  The  proving  ring  and  its  material  prop¬ 
erties  are  shown  in  Figure  7,  and  preliminary  results  with  PAFEC  level  5.2  are  shown 
in  Figure  8. 

The  initial  results  reveal  significant  discrepancies,  and  when  after  discussion  no 
satisfactory  reason  could  be  found  10  ,  a  detailed  examination  of  PAFEC  code  was  un¬ 
dertaken.  Various  differences  between  implemented  code  and  plasticity  equations  were 
found  (listed  in  Appendix  B),  most  of  them  having  only  minor  effects  on  results.  However 
changes  in  the  way  of  satisfying  convergence,  either  by  imposing  stricter  requirements, 
or  by  adopting  a  self-correcting  scheme  17  • 18  ,  produced  significant  changes  in  results 
for  this  problem.  It  is  expected  that  updated  levels  of  PAFEC  code  will  incorporate 
appropriate  modifications  16  . 

Alternative  FE  programs  such  as  NASTRAN  were  not  accessible  at  ARL  at  this 
time,  and  no  satisfactory  published  results  for  suitable  test  cases  were  available.  So 
in  order  to  test  the  modifications,  original  theoretical  solutions  were  developed  for  test 
problems,  given  in  Appendix  A.  Starting  with  the  simple  case  of  a  rectangular  beam  in 
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pure  bending,  solutions  are  extended  to  a  curved  rectangular  bar  in  pure  bending,  and 
finally  the  proving  ring  which  is  a  curved  bar  under  both  bending  and  tension. 

Results  for  the  straight  beam  and  curved  bar  in  pure  bending  are  shown  in  Figure 
9,  and  for  the  proving  ring  in  Figure  8.  In  all  cases  excellent  agreement  is  obtained 
between  the  modified  PAFEC  and  the  theoretical  solutions.  NASTRAN  had  by  this 
time  become  available  at  ARL,  so  it  was  used  for  the  proving  ring,  giving  a  solution 
virtually  coincident  with  modified  PAFEC. 

Cyclic  loadiug  for  the  proving  ring  is  shown  in  Figure  10,  comparing  isotropic, 
kinematic  and  mixed  hardening  (for  p  =  0.25),  with  results  qualitatively  similar  to 
those  given  in®  . 

6.  APPLICATION  AND  DISCUSSION 

Simple  seini-empirical  rules  for  dealing  with  plastic  deformation  at  stress  concen¬ 
trations  when  the  plastic  properties  of  the  material  are  known,  have  been  developed  by 
extending  solutions  obtained  for  particular  geometries  to  more  general  cases.  At  a  notch 
root  with  shear  stresses,  Neuber  23  derived  the  relation: 

K]  =  Ka  K,  or  K,  =  Kj  (24) 

where  A'/  =  theoretical  elastic  stress  concentration  factor 
K„  =  stress  factor  <t/<t 
A,  =  strain  factor  t /t 

Stowell 24  obtained  a  series  solution  for  elasto  plastic  deformation  around  a  hole 
in  an  infinite  plate  under  imiaxial  stress,  which  when  modified  to  cover  general  stress 
concentrations  25  is  expressed: 

K,  =  (K,-l)  Ko_  (25) 

Compressor  and  turbine  discs  subject  to  rotational,  other  mechanical  and  possibly 
also  thermal  loadiug  are  generally  the  most  critical  aircraft  engine  components  subject 
to  low  cycle  fatigue.  A  major  failure  source  is  at  holes  in  the  disc  web,  where  stress 
conditions  surroimding  the  hole  are  substantially  biaxial. 

For  a  circular  hole  in  an  infinite  thin  plate,  where  stress  conditions  resemble  those 
in  a  turbine  disc,  finite  element,  Neuber  and  modified  Stowell  solutions  are  compared. 
Principal  stress  ratios  rr2/rr t  =  0  (uniaxial),  0.5  and  1.0  (pure  biaxial)  are  shown  in 
Figures  11,  12  and  13  respectively.  Little  difference  between  the  three  stress  solutions 
can  be  seen,  largely  because  of  the  low  plastic  hardening  slope  (typical  Ti-8-1-1  material 
properties  were  used).  Strain  solutions  differ  significantly,  depending  on  the  principal 
stress  ratio.  For  uniaxial  stress  (<72  =  0,  Kt  =  3.0),  the  modified  Stowell  and  finite 
element  solutions  coincide,  at  least  for  lower  values  of  plastic  strain  (as  would  be  ex¬ 
pected,  since  the  Stowell  equation  derives  from  a  series  solution  for  this  case).  As  stress 
increases,  the  Stowell  equation  tends  to  overestimate  strain,  and  the  Neuber  solution  is 
similar,  excepting  that  it  overestimates  strain  throughout  the  stress  range. 

As  stresses  become  more  nearly  biaxial,  both  Neuber  and  modified  Stowell  rules 
continue  to  overestimate  strain.  With  pure  biaxial  stresses  ( rr2  =  cTi  ,  Kt  =  2.0). 
Neuber  and  Stowell  give  similar  moderate  overestimates  at  low  plastic  strains,  with 
Stowell  giving  much  higher  excess  estimates  of  strain  at  higher  stresses.  Fortunately 
such  estimates  lead  to  conservative  LCF  life  predictions  when  using  the  Coffin  Mausou 
law. 
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The  limitations  of  the  Neuber  and  modified  Stowell  rules  have  been  considered 
elsewhere21  .  Generally  speaking  they  overpredict  plastic  strain,  but  this  depends  on 
the  particular  geometry.  Most  investigations  have  been  concerned  with  their  accuracy 
for  LCF  fatigue  life  prediction 19  • 20  rather  than  for  strain  estimation,  and  in  this  respect 
they  are  widely  accepted  as  useful  design  rules. 

A  thin  disc  with  negligible  stresses  through  its  thickness  has  been  assumed  in  the 
current  analysis,  and  the  effect  of  biaxial  stresses  in  the  thickness  direction  at  the  hole 
surface  has  not  been  considered.  With  surface  biaxial  stresses  a  correction  should  be 
applied21  ,  which  tends  to  reduce  predicted  strains.  The  preliminary  applications  de¬ 
scribed  herein  are  intended  to  be  followed  by  a  fuller  explorations  of  thick  discs,  plastic 
crack  growth  and  other  practical  LCF  applications  in  the  next  phase  of  this  investigation. 

7.  CONCLUSIONS 

Application  of  PAFEC  to  cyclic  elasto-plastic  calculation  has  been  established  at 
ARL.  Careful  investigation  after  finding  discrepancies  between  results  and  published 
solutions  led  to  the  development  and  verification  of  modified  code.  This  code  extended 
the  available  strain  hardening  models  by  adding  mixed  hardening,  for  which  a  simple 
formulation  was  developed,  to  the  existing  isotropic  and  kinematic  hardening  rules. 
Verification  entailed  the  determination  of  original  alternative  solutions  to  problems  of 
bending  and  stretching  of  straight  and  curved  bars,  the  latter  exemplified  by  a  proving 
ring  for  which  several  alternative  FE  solutions  were  available. 

Sample  applications  were  made  to  holes  in  plates  with  biaxial  stress  similar  to  disc 
webs,  and  results  compared  with  the  Neuber  and  modified  Stowell  rules,  commonly  used 
for  design  life  estimation.  Compared  with  finite  element  solutions,  these  rules  tend  to 
overestimate  strain  in  biaxial  stress  conditions,  which  would  lead  to  a  conservative  life 
estimate. 

A  further  complexity  ensues  when  the  disc  or  plate  is  of  sufficient  thickness  that 
stress  levels  are  significant  in  the  thickness  direction.  Stress  conditions  at  the  hole 
surface  are  then  biaxial,  and  corrections  need  to  be  applied  to  give  equivalent  uniaxial 
strains.  Multiaxial  fatigue  is  a  complex  topic  and  is  not  considered  here;  it  has  been 
addressed  in  numerous  papers  elsewhere  26  .  Practical  discs  generally  exhibit  biaxial  and 
often  triaxial  effects,  and  analysis  of  these  discs  is  proposed  in  the  next  phase  of  this 
work. 
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APPENDIX  A  -  STRAIGHT  AND  CURVED  BARS 


Pure  Bending  of  Straight  Rectangular  Bar 

Adopting  the  standard  beam  assumption  that  plane  sections  remain  plane,  and 
using  notation  as  in  Figure  14,  strains  are  given  by: 

lx  g 

i  =  (mu  -r-  aad  I  =  — - 
*  ™  tmoz 

if  the  elastic  region  is  of  depth  g ,  yield  strain  ty  .  Stress  strain  relations  for  a  bilinear 
material  with  elastic  modulus  E ,  total  plastic  modulus  H  are: 

-ty  <  t  <  ty  a  =  E<  (26) 


(  >  £„  i 

v  a  =  (E  —  H)ty  +  Hi 

t  <  > 


Moment  M  per  unit  width  of  bar  is: 


M  =  2  /  crxdx 


Etmnj  ,  xdz  2 


/.'■«* 


H)iy  +  Htma/~}xdz 


At  initial  yield  My  =  E<y/i2/G  so  the  integral  becomes: 

M  1,.  .  H imaj  ,  , 

w  =,(1-  f)(3- --)  +  p  (28) 

My  2  h  E  ey 

Bending  and  Stretching  of  Curved  Bar 

Geometry  of  Deformation 

A  small  segment  d^  of  a  rectangular  bar  with  circular  centreline  is  shown  >n  Fig¬ 
ure  15,  and  the  following  geometrical  relations  apply: 

Total  rotation  of  face  w  =  da  +  id<f>  (29) 

where  c  is  uniform  strain,  da  is  rotation,  so: 

du)  da  , 

di  =  di+i  =  fi  +  i  (30) 

Rotation  of  the  bar  centreline  is  greater  by  the  increase  of  shear  deflection  from  7  to 
7  +  d7  ,  so  the  equation  for  change  of  curvature  is: 


so  since  i  is  small: 


r0{d^(l  +  i )  +  da  +  d7}  =  rd<f>(l  +  i) 


=  ia+^, 

r  r0  f(l  +  t)  W 


(31) 


Stresses  and  Strains 


Fibre  strain 


t  - 


rid<fr  +  (r  -  r)da 
rd<j> 


=  <+/*(  1-  f) 


(32) 


Yield  in  tension  and  compression  occnrs  where  t  =  ±  ty  ,  with  corresponding  radii: 

:i  - 


fir 

ft  +  t  ±  ty 


(33) 


Fibre  stresses  are  found  with  equations  (26)  and  (27): 

Elastic  region  r,  <  r  <  rt  rr  —  Et 


Plastic  region,  assuming  r/  >  r, : 

Tension  r  >  r( 
Compression  r  <  rt 


=  £{<  +  /?(!-  r)} 


<r  =  ±(£  -  H)t„  +  H{t  +  ft)  -  Hft 


(34) 


(35) 


Force  and  moment  per  unit  width  are  found  by  integrating  stress  over  the  cross  section, 
using  the  appropriate  equation  for  elastic  or  plastic  stress  regions. 

F  =  I  a  dr  (36) 


Jr  i 

=  f 

Jr, 


M  =  /  <t  rdr  (37) 

Explicit  integrals  for  radius  limits  «.  b,  can  be  found  from  expressions  (34)  or  (35): 

Elastic  Region  Forces  F  =  E\(ft  +  *  )r  -  ftr  log,  r]*  (38) 

Elastic  Region  Moments  M  =  E[(ft  +  t  )r2/ 2  -  /9rr j*  (39) 

Plastic  Region  Forces  F=  [^{  ±(E  -  H)ty  +  H(ft  +  i  )}r  -  Hftr  log,  r  (40) 

Plastic  Region  Moments  M  =  J{±(F  -  H)(y  +  H(ft  +  t)}r2/2  —  (41) 

Using  radius  limits  found  from  equation  (33),  values  of  F  and  M  for  given  ft  and  t  can 
be  found  from  equations  (38)  (41).  To  find  ft  and  i  as  functions  of  F  and  Af  requires 

solution  of  nonlinear  simultaneous  equations.  Since  the  diagonal  terms  predominate  (F 
is  principally  a  function  of  i  and  M  of  ft),  simple  linear  inverse  interpolation  can  be 
used.  (A  Newton-Raphson  solution  was  first  programmed  but  was  not  necessary.) 

Pure  Bending  of  Curved  Bar 

Shear  is  zero,  and  for  each  value  of  ft ,  a  value  of  <  is  found  giving  no  net  force 
(F  =  0) .  The  resultant  value  of  M  is  expressed  in  the  form: 

M/My  =  f(ft/fty)  (42) 

where  My  ,  fty  are  values  of  moment  and  rotation  at  initial  yield  (at  the  intrados). 


Circular  Proving  Ring 

For  symmetrical  loading,  the  ring  can  be  represented  by  one  quadrant  (Fig.  1C). 
Forces  and  moments  at  angle  <f>  for  load  2 P  are: 

Normal  Force  F  =  Pcos<f>  (43) 

Shear  Force  V  =  -Psin^  (44) 

Moment  M  =  Mo  +  Pr  cos  <)>  (45) 

Deflection  Equations 

(i)  Rotation 

Symmetry  requires  a  boundary  condition  such  that  there  is  no  net  rotation  over  the 
quadrant: 

/’  =  f  +  =  °  (4G) 

The  value  of  Mo  for  given  P  can  be  found  when  this  integral  is  satisfied. 

(ii)  Centreline  Deflection  of  Ring 

Geometry  of  deformed  ring  is  shown  in  Figure  17,  and  neglecting  second  order 
terms: 

.  .  .  ,  1  d<f> 

Initial  curvature  -  =  — 


Deflected  curvature 


1  _  d<j>  +  A  d<f>  _d<^  +J^ds 
n  ds  +  Ads  ds  (1  —  u/rj 


1  l  r  +  A 


1  d2u 

1  1  +  j"2  1 

Change  in  curvature  -  =  -- - — - 

rt  f  1  -  u/r  r 

Using  equation  (31),  and  noting  that  opposite  sense  of  M  in  Figures  15  and  16  reverses 
signs,  the  differential  equation  for  deflection  becomes: 

£+.-f(l -«)(#  +  £>  (481 

Shear  deflection  7  is  small  and  is  assumed  equal  to  the  elastic  deflection: 

kV 

1  ~  GA 

where  k  is  the  shape  coefficient  for  shear  and  A  cross-sectional  area.  Substituting  V 
from  equation  (44): 

_ kP  cos  <P 

d<t>  GA  '  ’ 

so  the  deflection  equation  becomes: 

d?u  wa  *Pcos<6.  .... 

w  +«-r0  -*)(/»-  -^  )  (SOI 

The  differential  equation  (50)  is  integrated  using  any  standard  numerical  integration 
technique,  with  terms  i  and  ft  for  each  angle  <f>  calculated  as  previously  indicated. 


APPENDIX  B  -  MODIFICATIONS  TO  PAFEC  CODE 


The  modifications  listed  were  applied  to  PAFEC  plasticity  code  at  level  5.2.  It 
is  expected17  that  later  releases  of  PAFEC  will  include  these  or  similar  modifications 
having  the  same  effect. 

(i)  Damping.  When  oscillatory  behaviour  occurs  and  convergence  is  uncer¬ 
tain,  a  damping  routine  may  be  invoked.  This  may  lead  to  incorrect  implementation 
of  the  Prandtl-Reuss  equations  and  has  been  discarded.  Section  (iii)  describes  other 
means  used  for  securing  convergence. 

(ii)  Kinematic  Hardening.  Implementation  of  the  Prandtl-Reuss  equations 
has  been  modified  to  properly  incorporate  yield  surface  translation. 

(iii)  Convergence.  Convergence  is  accelerated  and  stability  improved  by 
adaptive  adjustment  of  a  convergence  parameter. 

(a)  Convergence  is  accelerated  when  too  slow. 

(b)  When  maximum  error  increases  (instability),  previous  converged  values 
are  restored  and  convergence  slowed  until  error  reduces. 

(c)  After  convergence  is  re-established,  it  is  again  accelerated  to  hasten  at¬ 
tainment  of  solution  within  tolerance. 

(iv)  Iteration  Loop.  Exit  conditions  have  been  changed  so  that  a  self- 
correcting  load  adjustment 18  is  applied  at  the  start  of  the  next  load  increment.  This 
greatly  improves  solution  accuracy  for  given  error  tolerance,  and  may  allow  tolerance 
level  to  be  relaxed  while  maintaining  satisfactory  accuracy. 

(v)  Mixed  Hardening.  Mixed  kinematic  isotropic  hardening  with  variable 
proportioning  has  been  included. 
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FIG.  6  PAFEC  PLASTICITY  ALGORITHM 


FIG.  7  PROVING  RING  -  DIMENSIONS  AND  PROPERTIES 
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FIG.  8  COMPARISON  OF  SOLUTIONS  FOR  PROVING  RING  (SEE  FIG.  7) 
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FIG.  9  STRAIGHT  CANTILEVER  AND  CURVED  BAR  IN  PURE  BENDING 
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FIG.  13  HOLE  IN  PLATE  BIAXIAL  STRESSES 
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FIG.  14  RECTANGULAR  BAR  -  STRESSES  AND  STRAINS 
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FIG.  1 5  CURVED  BAR  -  DEFLECTION  AND  STRAIN 
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